transit_service_analyst demo

Import libraries

%matplotlib inline
import pandas as pd
import geopandas as gpd
import transit_service_analyst as tsa
import plotly.express as px
import shapely.geometry
import numpy as np
import plotly.offline as pyo
#import nbformat
#import plotly.io as pio

Create an instance of transit_service_analyst using a gtfs feed on file

path = r'R:/e2projects_two/Angela/transit_routes_2018/latest_combined'
transit_analyst = tsa.load_gtfs(path, '20180417')

c:\stefan\transit_service_analyst\transit_service_analyst\gtfs_service.py:32: DtypeWarning:

Columns (6) have mixed types.Specify dtype option on import or set low_memory=False.

The function get_lines_gdf() returns a GeodatFrame for each unique line by route, as determined by stop pattern. A trip_id representing each stop pattern is used as the unique ID and is called rep_trip_id. Show the first five rows:

transit_analyst.get_lines_gdf().head()

	shape_id	geometry	route_id	rep_trip_id	service_id	trip_headsign	direction_id	block_id	wheelchair_accessible	bikes_allowed	...	peak_flag	fare_id	agency_id	route_short_name	route_long_name	route_type	route_color	route_text_color	route_desc	route_url
0	CT_510:74	LINESTRING (-122.19550 47.97913, -122.19593 47...	CT_510	CT_7554558	1	Seattle 5th Ave	1.0	b_KPOB-CS:221:0:Weekday:1:18MAR:85001:12345	1.0	1.0	...	NaN	NaN	29	510	Everett - Seattle	3	0070c0	ffffff	NaN	NaN
1	CT_510:75	LINESTRING (-122.19773 47.97482, -122.19777 47...	CT_510	CT_7554559	1	Seattle 5th Ave	1.0	b_KPOB-CS:221:0:Weekday:1:18MAR:85001:12345	1.0	1.0	...	NaN	NaN	29	510	Everett - Seattle	3	0070c0	ffffff	NaN	NaN
2	CT_512:63	LINESTRING (-122.32897 47.59883, -122.32897 47...	CT_512	CT_7554560	1	Everett Station via Lynnwood	1.0	b_KPOB-CS:221:0:Weekday:1:18MAR:85001:12345	1.0	1.0	...	NaN	NaN	29	512	Everett - Seattle	3	0070c0	ffffff	NaN	NaN
3	CT_421:36:00	LINESTRING (-122.33432 47.61585, -122.33437 47...	CT_421	CT_7554561	1	Marysville	0.0	b_KPOB-CS:221:0:Weekday:1:18MAR:84050:12345	1.0	1.0	...	NaN	NaN	29	421	Marysville - Seattle	3	0070c0	ffffff	NaN	NaN
4	CT_511:73	LINESTRING (-122.32897 47.59883, -122.32897 47...	CT_511	CT_7554563	1	Ash Way P&R Lynnwood TC	0.0	b_KPOB-CS:221:0:Weekday:1:18MAR:85015:12345	1.0	1.0	...	NaN	NaN	29	511	Ash Way - Seattle	3	0070c0	ffffff	NaN	NaN

5 rows × 21 columns

Function to map line geometry

# Functiont to plot lines
def create_line_map(df, name_col):
    if len(df) == 0:
        print ('Empty DF. No lines to map!')
    else:
        lats = []
        lons = []
        names = []

        for feature, name in zip(df.geometry, df[name_col]):
            if isinstance(feature, shapely.geometry.linestring.LineString):
                linestrings = [feature]
            elif isinstance(feature, shapely.geometry.multilinestring.MultiLineString):
                linestrings = feature.geoms
            else:
                continue
            for linestring in linestrings:
                x, y = linestring.xy
                lats = np.append(lats, y)
                lons = np.append(lons, x)
                names = np.append(names, [name]*len(y))
                lats = np.append(lats, None)
                lons = np.append(lons, None)
                names = np.append(names, None)


        fig = px.line_mapbox(lat=lats, lon=lons, hover_name=names,
                     mapbox_style="open-street-map", zoom=8)
        fig.show()

Let's map one line for each route:

# Call function to map routes
geo_df = transit_analyst.get_lines_gdf()
create_line_map(geo_df.groupby('route_id').first().reset_index(), 'route_short_name')

The field transit_type is included so, for example, we can plot all routes that are not bus.

not_bus = geo_df[geo_df['route_type']!=3]
create_line_map(not_bus, 'route_short_name')

The function get_tph_by_line() returns a DataFrame with hourly frequencies for each rep_trip. Show the first 5 rows:

transit_analyst.get_tph_by_line().head()

	rep_trip_id	hour_2	hour_3	hour_4	hour_5	hour_6	hour_7	hour_8	hour_9	hour_10	...	hour_21	hour_22	hour_23	hour_24	hour_25	hour_26	hour_27	hour_28	route_id	direction_id
0	CT_7554558	0.0	0.0	3.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	CT_510	1.0
1	CT_7554559	0.0	0.0	0.0	4.0	6.0	4.0	4.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	CT_510	1.0
2	CT_7554560	0.0	0.0	0.0	3.0	4.0	4.0	4.0	4.0	4.0	...	3.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	CT_512	1.0
3	CT_7554561	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	CT_421	0.0
4	CT_7554563	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	CT_511	0.0

5 rows × 30 columns

Use get_tph_by_line() to get all rep_trips that have frequent service of 4 or more trips per hour during peak hours.

list_of_hours = ['hour_6', 'hour_7', 'hour_8', 'hour_15', 'hour_16', 'hour_17']
min_tph = 4
freq = transit_analyst.get_tph_by_line()
# create binary column for each hour
cols = []
for hour in list_of_hours:
    #print (len(freq))
    freq = freq[freq[hour]>=min_tph]

freq_rep_trips_gdf = geo_df[geo_df['rep_trip_id'].isin(freq['rep_trip_id'])]

print (f'There are {len(freq)} lines that have {min_tph} trips or more per hour.')

There are 53 lines that have 4 trips or more per hour.

# map them
create_line_map(freq_rep_trips_gdf, 'route_short_name')

The get_line_stops_gdf() function returns a GeoDataFrame consisting of the sequence of stops for each route returned by get_lines_gdf(). Use it and stops property to get the stops associated with the frequent routes.

# first get route stops for these lines:
route_stops = transit_analyst.get_line_stops_gdf()
route_stops = route_stops[route_stops['rep_trip_id'].isin(freq['rep_trip_id'])]

# Now get the unique stops:
unique_stops = route_stops[route_stops['stop_id'].isin(transit_analyst.stops['stop_id'])]


fig = px.scatter_mapbox(unique_stops,
                        lat=unique_stops['stop_lat'],
                        lon=unique_stops['stop_lon'],
                        hover_name="stop_id",
                        zoom=8,
                        mapbox_style="open-street-map")
fig.show()

The function get_tph_at_stops() function returns a DataFrame of total trips at each stop, independent of route. Let's use it to map stops by total number of trips for the AM peak:

list_of_hours = ['hour_6', 'hour_7', 'hour_8']
trips_at_stops = transit_analyst.get_tph_at_stops()
trips_at_stops['am_trips'] = trips_at_stops['hour_6'] + trips_at_stops['hour_7'] + trips_at_stops['hour_8']

# join to stops
trips_at_stops = trips_at_stops.merge(transit_analyst.stops, how = 'left', on = 'stop_id')

# get the stops with the most trips:
max_trips = trips_at_stops['am_trips'].max()
test = trips_at_stops[trips_at_stops['am_trips']==max_trips]

if len(test)==1:
    stop_id = test['stop_id'].values[0]
    number_of_trips = test['am_trips'].values[0]
    print (f'stop_id {str(stop_id)} has the most AM trips with {number_of_trips}')
    

stop_id KC_420 has the most AM trips with 263.0

# limit to 1000 stops
if len(trips_at_stops) > 1000:
    trips_at_stops.sort_values(by='am_trips', ascending=False, inplace = True)
    trips_at_stops = trips_at_stops.head(1000)

fig = px.scatter_mapbox(trips_at_stops,
                        lat=trips_at_stops['stop_lat'],
                        lon=trips_at_stops['stop_lon'],
                        hover_name="stop_id",
                        zoom=8,
                        mapbox_style="open-street-map",
                        size = 'am_trips')

fig.show()

Let's use GeoPandas buffering function to buffer these stops by 1/4 mile.

gdf = trips_at_stops.copy()
gdf = gpd.GeoDataFrame(gdf, geometry=gpd.points_from_xy(gdf.stop_lon, gdf.stop_lat))
# set initial geographic coordinate system.
gdf =gdf.set_crs(4326)
# need to switch to a projected system (state plane) to do buffer.
gdf = gdf.to_crs(2285)
gdf.geometry = gdf.geometry.buffer(1320)
# switch back to map
gdf = gdf.to_crs(4326)
gdf.set_index('stop_id', inplace = True)

fig = px.choropleth_mapbox(gdf,
                           geojson=gdf.geometry,
                           locations=gdf.index,
                           color="am_trips",
                           center={"lat": 47.6, "lon": -122.7073},
                           mapbox_style="open-street-map",
                           zoom=8)
fig.show()